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I. INTRODUCTION 


There have been many studies of geophysical fluid dynamics that deal with flow under 
the combined influence of buoyancy and rotation. The flows which are seen in cylindrical 
annulus experiments in which the thermal forcing and the rate of rotation are varied [1,2] range 
from simple axisymmetric flows to a full spectrum of waves interacting in a chaotic fashion. In a 
hemispherical shell space-laboratory experiment, Hart et al. [3] found a wide range of flow 
regimes depending upon the heating configuration and rotation rate. Such experiments are aimed 
at enhancing our understanding of atmospheric and oceanic dynamics by studying controlled, 
simple systems in which the flow dynamics can be isolated. However, any experiments are 
subject to instrumentation errors that affect the fidelity of the controlled conditions and to the 
data acquisition difficulties that prohibit detailed analyses of the phenomona under investigation, 
as well as to the expensive cost that discourages systematic and wider range parametric studies. 

A carefully developed numerical model that can simulate the experiments is a viable approach to 
the solutions of these problems. This approach allows more accurate control on the experimental 
conditions. It provides a complete data source for performing diagnostic studies. It can cost very 
little to run a high resolution model on a supercomputer if the model is coded efficiently. Our 
objective has been to develop an effective model that is able to simulate a broad range of 
configurations, including those mentioned above. 

The purpose of this paper is to describe the GEOSIM (geophysical fluid flow simulation) 
model which has recently been developed at the Marshall Space Flight Center (MSFC). It is a 
result of extensive renovation and expansion of the model used by Miller and Gall [4,5]. The 
model has been validated by comparing its results with those of laboratory experiments and other 
numerical models. Miller and Butler [6] and Miller and Fehribach [7] have given recent reports 
on the use of this model, although their papers use only part of the full capability of the model. 

In this paper, we present other validation results which are of interest to the community. These 
cases are: (1) the steady baroclinic wave in a rotating side-heated annulus [8], (2) the amplitude 
vacillation studies of Buzyna et al. [9], (3) transition to baroclinic waves in the bottom-heated 
annulus.experiments at MSFC [10], and (4) the Spacelab Geophysical Fluid Flow Cell (GFFC) 
observations described by Hart et al. [3]. 


n. DESCRIPTION OF THE MODEL 
A. Equations 

The model is based upon the Boussinesq Navier-Stokes equations in rotating spherical 
coordinates. The equations, with the advective terms written in flux form, are as follows. The 
definitions of the symbols are given in Table 1. 





Table 1. List of Symbols. 


Inner radius 
Depth 

Body force ("gravity") 

Gravity at r = a 
Longitudinal wavenumber 
Pressure 

Radial coordinate, distance from center 
Time 

Temperature 

Longitudinal velocity (positive for westerly flow) 

Meridional velocity (positive for southerly flow) 

Radial velocity 

Thermal expansivity 

Latitudinal coordinate 

Thermal diffusivity 

Longitudinal coordinate 

P/ p 0 (often referred to a "pressure" in text) 

Reference density 

Kinematic viscosity 

Rotation rate 


2 




( 1 ) 


du ( 1 3 , xx . 1 3 , J2, . 1 du 2 \ . uv tan <|> «w , A 

■3T =- (775F* 5$ <“ v C0J ♦> + 7 3F (M, “^ )+ 77^* STJ + ~ r ~ T +2ftv *" ♦ 

- 2£lw cos <|> + v[v 2 h - - (« + 2^ (v «h #- w «w ♦))] - ^^ 


3v / 1 3 / 7 x \ . 1 3 / _ _7 x . 1 3uv\ h tan <|) vw 

W = " ^ ^ 3F {wv * )+ TESTS ~ST) r T 

. , r o2 v 2 5 wi 6 3 u , 2 3 wl 1 3 n 

-2fta «n 4» + v|^V V - ^ + JT^T^ 3X + ^ ^ J ~ 7 W 


( 2 ) 


3w ( \ 3 , xx . 1 3 / 2 ji\ . 1 3hw\ , (h 2 + v 2 ) 

W = ~ (TESTS ss ^ C0S ^ + TTr ^^ + 7ESTSSr) + ~ ~ 

, fro 2 3 , .. 2 du 2 wl , „ _ 311 

+ 2 fl« cos S + v [ v 2 W - jr^ ^ (v a» +) - ^-^rj +P^- -gp 


( 3 ) 


§ = ~ (FE5F? |f (vr « + 7 I (w7>a) + FSF? THf) + K(v2r) * (4) 

The continuity equation is as follows:. 
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The spherical three-dimensional Laplacian is: 
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B. Numerical Method 

The model uses a mixed spectral and finite difference technique, similar to that used by 
Quon [11] for a cylindrical model. In the azimuthal direction, each field is decomposed into 
Fourier components of the form A t exp(/£X,), where "K is longitude, k is a given wavenumber (0 for 

axisymmetric), and A k is the complex amplitude of one of the dependent variables. By 
substituting the Fourier series for each variable into equations (l)-(4) and operating on both sides 
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by the orthogonalization procedure, one obtains a set of equations for each Fourier mode [11]. 
These equations can be simplified in a consistent manner to represent various mathematical 
systems in different degrees of complexity. GEOSIM is able to solve five classes of problems, 
which are handled consistently with the same numerical algorithms: 

1. Axisymmetric flow 

2. Linear instability problem with respect to a previously computed, fixed axisymmetric 

state 

3. Linear instability problem with respect to an evolving axisymmetric state as in 1. 

4. Nonlinear problem for single and multiple waves, in which the only quadratic terms 
retained are those due to the interaction between the wave and axisymmetric parts of the flow 
(hereafter called the wave-mean flow interaction case) 

5. Fully nonlinear problem for multiple waves with complete wave-mean and wave-wave 
interactions. 

Note that modes 2 and 3 compute wave structure and its corresponding growth rate which 
asymptotically becomes the fastest among the modes of the linear system. The detail of 
computing the eigenfunction and the eigenvalue is given in the next paragraph. Mode 4 allows 
an arbitrary sequence of wavenumbers in the model. Mode 5 deals with a simple spectral 
truncation series with wavenumber being multiplied by a fixed harmonic factor. Therefore, in 
addition to the traditional spectral resolution 1 to K, the model can solve for the same number of 
wave components but in a higher harmonic, e.g., wavenumbers 2, 4, ..., 2K. In mode 5 
application, the well-known transform method (see ref. [12]) is employed to compute the spectral 
transformations of nonlinear terms in the equations. Details of this method are given in section 
ILF. 


The method of computing the largest real eigenvalue and its associated eigenfunction 
deserves special attention. The model is integrated forward in time in the expectation that the 
solution will eventually become dominated by the fastest-growing (or slowest-decaying) 
eigenmode. At each time step an estimate of the linear, complex growth rate (co) is calculated for 
the wave by spatially integrating the temperature tendency terms. The term (-0)0 is added to 
each of the wave prognostic equations (where Q refers to the prognostic variable). If the fields 
become dominated by a single fastest-growing (or slowest-decaying) eigenmode, this additional 
term stops the wave from growing (or decaying) and propagating. For most situations, this 
indeed occurs, and the system of equations (with the additional term -(oQ) has a steady state 
solution. The use of this technique avoids the situation of the model "blowing up" due to 
exponential growth, it makes the solution independent of the time-stepping formulation, and it 
allows the use of differential time steps to enhance the convergence rate. 
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A similar procedure can be applied as an option to obtain a steady state solution for the 
mode 4 application when the solution includes a wave of constant structure which propagates 
azimuthally. In that case, only the imaginary part of the growth rate is used in the additional 
term. After a steady wave is achieved in the solution, this procedure stops the wave from 
propagating in the azimuthal direction and permits the calculation of a time-independent solution. 

In the cases in options (1) through (4) where the steady state is the desired result, the use 
of differential time steps (different values of the time step for different equations) can 
dramatically speed up the rate of convergence. In particular, this is true for cases in which either 
momentum or heat (usually the latter) has a much lower diffusivity than the other. The time step 
constraint in the cases of interest is often due to the diffusive terms; hence, the time step for the 
more slowly diffusing component can be made several times that of the other component without 
loss of numerical stability or accuracy (for steady state flows). Of course, variable time steps are 
never used in cases where the time evolution is important (for example, vacillation studies). 

C. Spatial Discretization 

1. The Grid System . The model uses a staggered grid system (see Figure 1) similar to 
that of Quon [11]. The physical domain is partitioned into an arrangement of J by K (JJ-2 by 
KK-2, in the code’s nomenclature) elements in the meridional and vertical directions. The 
variables u, T, and n are defined in the center of the elements; v is defined at the middle of the 
vertical boundary of each element (i.e., on the sides); and w is defined at the middle of the top 
and bottom of each element. Note that the grid system includes artificial points outside the 
physical domain, for efficiency in computation and accuracy in imposing the boundary 
conditions. The grid can be stretched in both directions, so that finer resolution can be obtained 
in regions of strong gradients. The initialization routine of the model constructs a stretched grid 
which has grid spacings that are determined based on the formula sin(jry) + 5, where y ranges 
from zero to one (proportional to physical distance) and 5 is an input parameter. For larger 8 the 
grid stretching is less severe. There is also an input parameter which adjusts the location where 
the grid is to be stretched. This is done by moving the sine function so that the desired amount of 
stretching in the location of interest is achieved. After the spacing function is determined, the 
grid intervals are multiplied by the appropriate constant such that they are in the proper units. 

The method differs from that of Quon, who transforms the equations to a stretched coordinate 
system and then discretizes with a constant interval. The advantage of the current method is to 
achieve increased flexibility, which is important because of the requirement to use the same 
model for a variety of geometries and flow types. 

2. Finite Differencing. The spatial derivatives in the Fourier transformed equations are 
approximated by a control volume discretization procedure. Standard centered difference 
approximations are used for the derivatives in each of the equations with the following 
exceptions. 

For the temperature equation, the terms representing advection of heat by the mean 
meridional flow are evaluated as a weighted average of upwind and centered differences. The 




Figure 1. Staggered grid system, showing the placement of the velocity components. The T and 
TI variables are placed on the same points as the u component. 
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relative weights given to the two differences is input by the user. A heavier weight on the 
upwind scheme may help suppress oscillatory behavior of the numerical solution due to the 
splitting property inherent in the centered differencing scheme. However, too much upwind 
differencing introduces artificial viscosity (see ref. [13]), which may degrade the accuracy of the 
solution. Thus, the requirements for upwind differencing must be evaluated for each individual 
problem. In general, we have used upwind differencing as little as possible, introducing the 
smallest fraction necessary to avoid small-scale features which are clearly numerical (which 
occur due to the advective time scale becoming smaller than the diffusive time scale). 

The second exception to the standard centered differencing formulation deals with the 
evaluation of the flux terms in the v and w equations. In that case, the Wam-Vamas et al. [14] 
"Scheme B" is used in order to conserve kinetic energy while retaining second order accuracy. 

D. Time Discretization 

A two-level scheme is used for the time differencing, which minimizes the computer 
memory requirements. First-order forward differencing is generally used, although exceptions to 
this rule are made according to numerical stability criteria. Those deviations are described in the 
following paragraphs. 

In order to lessen the time step constraint due to the presense of inertial gravity waves, the 
Coriolis and the body force terms in the w and v equations use newly updated values of u and T 
(i.e., u and T are updated first, before the tendencies of v and w are calculated). While this 
introduces no complication for the axisymmetric part of the flow where u is not involved in the 
continuity equation, mass continuity requires that all three velocity components be updated 
simultaneously in the wave equations. In the latter case, a tentative value of u, based upon the 
pressure at the previous time step, is used in the Coriolis terms for v and w. Experience has 
shown that the difference in the results obtained between this update procedure and the standard 
forward differencing scheme is insignificant, while the gain in allowable time can be sizable 
when the rotation rate or the vertical stratification is large. 

To further permit longer time steps, there is an ADI (alternating direction implicit) option 
for the diffusion terms and those involving advection by the mean flow. The ADI method treats 
terms implicitly in one direction and explicitly in the other direction for even time steps, and on 
the odd time steps switches directions. The following is a description of the ADI method applied 
to the vertical direction at odd time steps, using the T equation as an example. Letting n + 1 and n 
be the forecast and current time levels, respectively, the equation is discretized as follows: 

w " 5 ' = ~ ~ = ~( Dr 7 ’ ,,+1+ D * r " + Dx r ") + K ( L ’ r " +1 + L \ r " + *1 r " +1 )> < 7 ) 

where D^, D^, and D r represent the centered finite difference approximation of flux divergences 

due the axisymmetric components of u, v, and w, respectively, and the L 2 ' s represent the finite 
difference analog to the Laplacian components. The velocity components (valid at time level n) 
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are included in our definition of the operator D. The equations are rewritten in terms of AT and 
7*. This is done by adding and subtracting the terms on the right-hand side that are at time step 
n, then replacing (7" +1 - 7") with AT and gathering all the terms involving AT on the left-hand 
side: 

[ 1 - A/ (k (l) +lJ)) - £> r ] AT = discretized form of the (8) 

right-hand side of (4) 
involving values at time step n. 


When (8) is expanded into matrix form, it is a tridiagonal matrix system for AT and is 
solved by Gaussian elimination. At even time steps, the r-derivatives on the left-hand side are 
replaced by the (^-derivatives. Note that the X part of the diffusion is implicit and the X part of 
the advection is explicit at all time steps. 

E. Poisson Equation for Pressure 

At each time step, a pressure field must be found which results in mass continuity (zero 
divergence) for the new velocity field. The procedure is similar to that of Williams [15]. A 
Poisson equation is derived by first writing equations (l)-(3) in vector form: 

^ =- vn + c, (9) 

where V is the velocity vector and G represents the terms on the right-hand sides of equations 
(l)-(3) other than the pressure gradient terms. 

Replacing the left-hand side of (9) by forward differences and multiplying both sides by 
At gives: 


V* +1 - V" = -(Vn* + G*) At. (10) 

n* is the change in pressure from the previous time step and the gradient of the old n is included 
in the definition of G*. Applying the divergence operator on (10), setting the divergence of the 
new velocity vector to zero, and rearranging terms gives: 

(ID 

Note that the term V • V is kept on the right-hand side since there is some round-off divergence 
from the previous time step. The boundary condition used is that the normal derivative of n is 
zero. This is consistent with our formulation of the Laplacian, in which the normal derivative 
across the boundary is taken to be zero (justifiable because the momentum equation is not used to 
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predict the boundary values). As explained by Williams [15], this condition results in IT being 
actually a quantity equal to the physical pressure correction except at the boundaries. A 
correction to n* could then be applied in order to recover the physical quantity; however, that 
step is not necessary, since the boundary derivative of II is not required in solving the 
momentum equations. The elliptic equation with the Neumann boundary condition is solved 
using a direct method that is vectorized for the CRAY-XMP [16]. The solver has no problem 
with the fact that there is not a unique solution for the axisymmetric part. It should be noted that 
the solution of the pressure equation is the most time consuming part of the model. For this 
reason, we provide a control that skips the pressure update if the maximum absolute value of the 
divergence is smaller than a prescribed value. 

F. FFT - Method for Solving for Nonlinear Terms 

Around an azimuthal ring, every dependent variable is expanded in a discrete complex 
Fourier series as follows: 


K 

q(kj)= Z Q extfkj), 

*=>-A * 


(12) 


where A .■ = InjlL, L is the number of intervals in the azimuthal direction, and K is the truncation 
limit of wavenumber k. The Fourier (or spectral) coefficients, Q k (r,§,t), are computed by the 
following: 


, t 

Q k = T = Zq(k)Qxp(ik\j). (13) 

^ j=\ j 

Note that (13) is exact only if the product q(kpexp(ikXj ) is a truncated.trigonometric series with 
maximum wavenumber less than or equal to L-l. Otherwise, aliases would result. Hence, if 
q(kp represents a product of two real variables representable by the same spectral series with 
truncation ‘K, then L > 3AT+1 must be satisfied. The actual calculations of (12) and (13) are 
carried out by a vectorized fast Fourier transform (FFT) algorithm [17]. This FFT is the same as 
that used by the NCAR Community Climate Model which limits L to have multiplication factors 
of 2, 3, and 5 only. This does not place a constraint upon K, although since optimal values for L 
as a function of the number of waves are tabulated in the present program, the number of waves 
must be <40. This number can easily be increased by expanding the table. 

Since q(kp is real, the complex discrete Fourier transform is conjugate symmetric; i.e., 

= Q*^. Hence, only k = 0 to K Fourier coefficients are calculated. In fact, (12) and (13) can 
be rewritten in a real trigonometric form as follows: 
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where 


2 r 

a* = ^Z q(kj) cos k\j 

(15) 

2 L 

bk = £ Z q(kj) sin kXj. 

(16) 


Hence, for each complex coefficient, Q# the real and imaginary parts are one half of the 

coefficients of cosine and sine series, except for the k = 0 component. The Fourier transform 
method described above (see also ref. [12]) is efficient in dealing with the nonlinear (quadratic) 
terms in the governing equations (l)-(4). In short, all variables required in the calculations are 
first synthesized into physical space, where the quadratic factors (e.g., uv, v 2 , vT, etc.) in the 
nonlinear terms are calculated. The spectral coefficients of those terms are calculated using the 
FFT. These spectral coefficients are then used in the finite differencing scheme to obtain the 
nonlinear tendencies for each wave component. Note that the FFT is used only when a fully 
nonlinear run of the model is made (mode 5). Otherwise, the nonlinear terms (which do not 
involve wave- wave interactions) are calculated in spectral space, on the r-< j) grid. 

G. Boundary Conditions 

The model fluid is confined in a latitudinal channel of arbitrary width (which may include 
the pole as part of the fluid domain) and depth around the polar axis. The use of latitudinal walls 
at any latitude is an advantage of the hybrid scheme; this would not be possible if, for instance, 
spherical harmonics were used. Each of the four computational boundaries has the option of 
being free-slip (zero stress) or no-slip (zero velocity), and each boundary may have a heated 
(fixed temperature) or adiabatic (zero normal gradient) temperature boundary condition. The 
temperature profile imposed upon a differentially heated horizontal boundary can have a linear, 
sine, or logarithmic shape, where the latter is defined by Equation (7) in Miller and Fehribach 
[7]. The vertical boundaries for the temperature have a linear profile. It is straightforward for a 
user to program a different temperature profile. The values of vertical velocity (w) on the top 
and bottom are set to zero regardless of the boundary condition configuration. The value of the 
v-component is set to zero on the two meridional boundaries, except for the case in which the 
domain extends to the pole (next paragraph). 

Polar Cases . When there is a pole in the domain (e.g., a full cylinder or spherical shell), 
the computional boundary condition there is non-trivial. The pole is a singular point in the 
coordinate system (i.e., the velocity components are not well defined there), but there are no 
physical constraints on the flow. Fortunately, the grid is structured such that a polar boundary 
condition is required only for the v-component of velocity. The values for the other variables are 
arbitrary because they are always multiplied by zero in the finite difference scheme. This is due 
to the fact that the differences are in flux form and the area of the pole itself is zero. However, 
the polar value of v is used in the calculation of viscous stress one-half grid interval from the pole 
in order to calculate the viscous tendency for v at the point one grid interval from the pole. By 
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considering a small neighborhood of the pole, it may be easily shown that the only Fourier 
component which includes non-zero flow at the pole is the wavenumber one component. 
Therefore, a correct condition for the axisymmetric and all wave components except 
wavenumber one is to set v zero at the pole. For the wavenumber one component, the correct 
boundary condition is zero meridional gradient. This arises from the fact that the flow resulting 
from the wavenumber one component along a great circle passing through the pole must be 
symmetric about the pole and must have a finite (Cartesian) derivative at the pole. This 
derivative is approximated by setting the polar value equal to that one grid point from the pole, 
which is only first-order accurate. 

To produce smooth cross-section plots, polar values of the other components are needed 
when plotting. A zero meridional gradient is used for the axisymmetric components of w and T, 
and zero amplitude values are used for the wave components. The boundary condition on u is 
taken to be the same as for v. 


fll. MODEL VALIDATION 
A. Axisymmetric Flow, Linear Waves, and Hysteresis 

Verification of the model was performed in the order of complexity of the problem. 
Axisymmetric results were compared and agree well with previous results by Williams [8], Quon 
[11], and Miller and Gall [4,5]. (Note that, while the computer code was rewritten, there was no 
fundamental change in the algorithm for the axisymmetric case from the Miller and Gall papers ). 
Linear wave results were compared with Miller and Gall [5]. Exact agreement in the latter case 
was not obtained nor expected, since the Miller and Gall model made the hydrostatic assumption 
in the wave component of the flow. The next step was to verify the model with the addition of 
the nonlinear wave-mean flow terms. Miller and Butler [6] used the model in mode 4 to 
calculate hysteresis in the transition curve as found experimentally by Fein [1], and good 
agreement with the experiment was obtained, as shown in Figure 2. 

B. Steady Baroclinic Annulus Wave: Comparison with Previous Work 

Attempts were made to validate a steady wave case with Williams [8] and Quon [11]. 
Note that Quon’s model is essentially the same as our model in the single wave-mean flow 
interaction mode, while William’s model is similar to our model in a multiple wave, fully 
nonlinear mode. With the periodicity condition applied to one-fifth of the annulus, the wave 
components resolvable in the Williams model are 5, 10, 15, etc. Figure 3 shows the azimuthal 
mean state from the GEOSIM model in mode 4, including only wavenumber 5. In comparison 
with the Quon and Williams results (which share resemblance), the wave has smaller amplitude 
and the mean fields are not affected as much by feedback from the wave. However, the 
differences reduced considerably when the model was run in fully nonlinear mode with waves 5 
and 10 included. When wave 15 was also included (Figure 4), we obtained results that virtually 
agreed with those of Quon and Williams. We also ran GEOSIM in single-wave mode with wave 
4 instead of wave 5 selected, and found that the resultant axisymmetric state and wave amplitude 
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OMEGA (sec- 1 ) 


Figure 2. Amplitude of the temperature wave as a function of rotation rate as determined by Fein 
[1] (dashed curve) and Miller and Butler [6], As the rotation is increased with an axisymmetric 
state established, the model predicts transition to baroclinic wave flow at point A. As the 
rotation rate is decreased with a wave established, the model predicts transition to axisymmetric 
flow at point F. 
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L - 


Mean State 

Zonal Velocity Maximum * 0.476 

Contour Interval » 0.05 Minimum * -0.082 



Mean State 

Stream Funoon Maximum * 0.038 

Contour Interval =* 0.005 Minimum * -0.0007 


Figure 3. Temperature, zonal velocity, and meridional stream function from the azimuthal mean 
state computed by GEOSIM in mode 4 for the case of Williams [8] when a single wave (5) is 
used. 
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Temperature 
Contour Interval » 0,5 


Mean Slate 
Maximum = 5.0 
Minimum =* 0.0 


Zonal Velocity 
Contour Interval = 0.05 


Mean State 
Maximum * 0.291 
Minimum =-0.105 



Figure 4. Similar to Figure 3, but computed by GEOSIM in mode 5 (fully non-linear) with 
waves 5, 10, and 15 included. 





were similar to the previous works. The discrepancy between our single wave case and that of 
Quon is left unexplained. 


C. Vacillation 

Wave amplitude vacillation is characterized by a periodic growth and decay of a 
baroclinic wave with a similar periodic decrease and increase of azimuthal mean baroclinicity. 
GEOSIM in mode 5, with wavenumbers 1 through 15, was run using the same geometric and 
dynamic parameters as one of the experiments of Buzyna et al. [9]. Figure 5 shows the time 
variations of kinetic energy integrals for the axisymmetric and wave components in about two 
vacillation cycles after 21,400 seconds (or approximately 600 rotations) of model integration 
time. As observed in the laboratory, wave 4 dominates after the system equilibrates to a steady 
vacillation. Wavenumber 8 and higher harmonics are present, but there is no side-band 
contribution to the energy spectrum. The period of vacillation is estimated to be 45 rotations, 
which is two-thirds of what was seen in the laboratory experiments. Figure 6 shows the 
temperature distributions of the symmetric state and wavenumber four disturbances at the middle 
level of the model fluid. The four diagrams represent the synoptics at the four phases of a wave 
amplitude vacillation cycle, denoted by WL, SH, WH, and SL. About 0.3 °C difference in 
temperature variation between a lowest wave (WL) and a highest wave (WH) stage can be 
detected, in reasonable agreement with the observations. The movement of the arrow mark 
indicates that waves travel around the annulus at a frequency about two times as fast as the wave 
amplitude changes, consistent with the experiments. A more detailed analysis of this and other 
numerical experiments will be presented in a future paper by Lu, Miller, and Butler. 

D. Bottom-Heated Cylindrical Annulus 

Laboratory experiments in which a temperature gradient is imposed on the lower surface 
of a cylindrical annulus have been performed at MSFC. The outer vertical wall of the annulus is 
also constrained to have a fixed temperature profile, while the upper and inner walls are made of 
a relatively non-conducting material (Plexiglass). Many of the results of these experiments are 
similar to those with the more traditional, side-heated annulus. For example, a boundary 
representing the transition between stable, axisymmetric flow and that with baroclinic instability 
in thermal Rossby number-Taylor number parameter space has a similar shape to the side-heated 
and -cooled experiments. The model was used to predict the transition curve, including number 
of the initial waves. Figure 7 shows the comparison between the laboratory experiments and the 
calculated transition. The agreement in placement of the transition curve is fair, but it is not as 
good as a similar comparison in the side-heated and -cooled case. The agreement in terms of 
wavenumbers is excellent. The authors believe that the differences in the placement of the 
transition curve is due to imperfect thermal boundary conditions in the experiment. Significant 
improvement from a previous comparison was obtained by adding more insulation to the 
apparatus in the laboratory experiments. 
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Figure 5. Time variations of kinetic energy integrals for the axisymmetric and wave components 
as calculated by GEOSIM for a vacillating case of Buzyna et al. [9], see case B-j of Pfeffer et al. 
[19]. 
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Thermal Rossby Number 



Figure 7. Calculated and experimental transitions for the bottom-heated annulus of Miller and 
Reynolds [10]. 
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E. Spherical Shell 


Comparisons have been made of fully three-dimensional flow patterns predicted by the 
model with experiments performed on the Spacelab 3 mission with the Geophysical Fluid Flow 
Cell (GFFC). Agreement was excellent to the extent that the experimental data allow 
comparison (which is limited to flow patterns). Documentation of these comparisons will appear 
in a separate paper by Miller and Leslie. More detailed comparisons can be made between the 
present model and another one, namely a model based on spherical harmonics developed by Gary 
Glatzmaier of Los Alamos National Laboratory [18]. Note that the Los Alamos model is 
constrained to have a free-slip condition at the equatorial wall, so the present model was run with 
that condition as well. An initial condition was determined by integrating the Glatzmaier model 
for a short time (see Figures 8 and 9). The fields at that point were then used as initial conditions 
for both models. After 9 minutes, the results from the two models are compared in Figures 8 and 
9 for temperature and for radial velocity at mid-depth. Note that the agreement is remarkable, in 
spite of the fact that the patterns have changed significantly from the initial conditions. Further 
comparisons are in progress for a case with stronger gravity, in which the fields change more 
quickly and in which the results have been shown to be extremely sensitive to initial conditions. 
These comparisons will be published separately by Johnston and Miller. 

IV. SUMMARY 

The MSFC-developed model GEOSIM has been shown to be an extremely powerful tool 
for simulating a wide range of geophysical flow problems. The model employs finite differences 
in the meridional plane and spectral decomposition in the azimuthal direction, which results in a 
code which is flexible in terms of domains and boundary conditions. The use of a spectral 
transformation in the azimuthal direction allows the calculation of axisymmetric states, linear 
waves, and states with wave-mean flow interaction only. The code has been validated by 
comparing its results to both experimental and numerical studies. The example cases that are 
given in this paper demonstrate the types of fluid flow problems which GEOSIM was designed 
to solve. Further research at MSFC using the model is planned, and the authors welcome 
inquiries from researchers with similar needs and interests who may wish to use the model. 
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TEMPERATURE AT 2.78 CM MAXIMUM = 10.17 

CONTOUR INTERVAL = 2.0 MINIMUM = -3.90 



TEMPERATURE AT 2.78 CM MAXIMUM = 11.08 

CONTOUR INTERVAL = 2.0 MINIMUM = -5.37 



TEMPERATURE AT 2.78 CM MAXIMUM = 11.03 

CONTOUR INTERVAL s 2.0 MINIMUM a -5.34 

Figure 8. Temperature fields from the comparison between GEOSIM and the model of 
Glatzmaier [18]. The upper plot is that of the initial condition (same for both models). The 
middle plot is from the GEOSIM model 9 minutes later, and the lower plot is from the 
Glatzmaier model at the same time. 


20 



VERTICAL VELOCITY AT 2.78 CM MAXIMUM =0.0087 

CONTOUR INTERVAL = 0.002 MINIMUM =-0.0082 


VERTICAL VELOCITY AT 2.78 CM MAXIMUM = 0.016 

CONTOUR INTERVAL = 0.005 MINIMUM = -0.014 




VERTICAL VELOCITY AT 2.78 CM MAXIMUM = 0.01 6 

CONTOUR INTERVAL = 0.005 MINIMUM = -0.015 

Figure 9. Similar to Figure 8, but for radial velocity component. 
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APPENDIX A: User’s Guide 


The current version of the GEOSIM model (1990) was developed and runs on the CRAY 
XMP at Marshall Space Flight Center (MSFC). While this User’s Guide often specifically refers 
to the Cray, we have generally avoided using instructions that are unique to this machine. With 
the construction of proper operating system commands, the model should run at any site that 
supports a FORTRAN 77 compiler. For example, it has been successfully run on an IBM RISC 
6000. 


Compilation Procedures 

The model’s source code is contained in four ASCII files: one for the main program and 
subroutines, one for subroutines used in the FFT and the elliptic equation solver, one for the 
derived parameter statements and common blocks, and one for the key parameter statements. 

The file with the common blocks is "included" in the main program and most of the subroutines 
at the time of compilation. The sizes of the work arrays are set by the user to the desired value in 
the key parameter statements. In the current implementation on the Cray, these parameter 
statements are in the same file as the operating system commands that compile the code. Note 
that, while the radial and latitudinal grids used must be the same as the array dimensions, the 
number of waves can be less than or equal to the compiled dimension. The following key 
parameters must be set at compilation: 

MIIW: Maximum number of wave components (excluding the axisymmetric 

component) 

JJ: Number of grid points in the meridional direction. Note that the number 

of interior grid intervals in JJ - 2. 

KK: Number of grid points in the vertical direction. The number of interior 

grid intervals is KK - 2. 

NONCASE: Selection for mode of operation (determines whether FFT is used). 

= 1 : Fully nonlinear run 
= 0: Any other mode of operation 

IWDIM: Dimension of the work array for the elliptic equation solver. A formula 

is given for estimating the minimum value of this parameter. If the 
value is too small, an error message is printed and the model does not 
run. 


Restart Files and Restart Procedure 

At the end of each integration, each wave component (excluding the axisymmetric part) is 
written to a separate logical unit (7 and above), and all waves are concatenated and written to a 
single logical unit (99). (Only one of these sets of files— whichever is more convenient for the 
user— needs to be saved.) The axisymmetric data contain grid information, are slightly different 
in format from the wave data, and are written to logical unit 1. The user can choose to save 
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whichever files desired for analysis and/or restarting the model. There axe several ways to 
initialize the model: 

ISTART = -2: Restart from scratch - This uses boundary temperature and 
profile information within the source code to set up the initial axisymmetric 
temperature field. (The interior temperature field can either be a constant value 
or an interpolation of the boundary values. This is determined by the input 
parameter "TINIT," which is on the last line of the input deck.) If waves are 
included, they are initialized with a small temperature disturbance (randomly 
oriented in the azimuthal direction) located in the middle of the meridional 

grid.* 

ISTART = - 1 : Axisymmetric restart - A previously computed axisymmetric 
field on unit 1 (with grid information) is used as an initial condition. The 
waves are initialized in the same way as a "scratch" run. 

ISTART = 1 : Both axisymmetric and wave fields are read from units 1 and 7 
(and up, if necessary), respectively. Note that the wave fields must be in 
ascending order 


*There is a call to the random number generating routine RANF in the subroutine INITIA. This 
is a spot where the user may need to make a change for a different computer installation. 
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Input Parameter File 


The unit 5 input file (usually concatenated with the operating system commands to run 
the model) prescribes various numerical and physical parameters necessary to run the model. 
Table A. 1 defines each of the input parameters. The grid and boundary temperature parameters 
are used only in the "scratch" starts; other parameters are used in all run modes and replace any 
values stored in the axisymmetric save file. For example, if an axisymmetric field is run with 
omega = 2.0 and restarted with omega = 1.8 in the runstream, the value of 1.8 is used in the new 
run. Unless a "delta omega" of -0.2 is input, the model does not add the differential spin-up 
rotation; that is, the u-field in the input file is taken to be relative to the new rotation rate. 
Boundary conditions on temperature cannot be changed in a restart. (It would be straightforward 
for the user to write a program which edits the restart file to change them, if desired.) 

Diagnostics such as integrated kinetic energy are calculated at intervals.which depend on 
the "NSTEP" and "ISKIP" factors. If the number of steps (NSTEP) is set to 5000 and the skip 
factor (ISKIP) is set to 1000, diagnostics.are be written to the output file on steps 
0,1, 1000,2000..., 5000. If ISKIP is set to zero, the default factor of NSTEP/200 is used for 
ISKIP, and therefore the diagnostics are written 202 times. A summary page, which contains 
basic information for the axisymmetric state and all of the waves, is written after every 5 pages 
of diagnostics. 

The "IRUNCD" option corresponds to the mode of integration that is to be run in the 
model. A "1" indicates an axisymmetric only run, "2" - linear waves only, "3" - axisymmetric 
and linear waves (no interaction), ”4” - axisymmetric and nonlinear wave (wave-mean flow 
interaction), and "5" - fully nonlinear wave run. The value of IRUNCD affects the interpretation 
of the "WAVE-NUMBERS" line. In modes 2 through 4, this refers to the wavenumbers 
included in the calculations. In mode 5, this line lists the waves to be included in the printed 
diagnostic output, the "WAVE FACTOR" indicating which set of wave numbers to run. For 
example in mode 5, if the number of waves is set to 10 and the wave factor is set to 3, then waves 
3,6, 9... 30 would be included in the calculation. Table A.2 gives a list of memory and time 
requirements for each ran option for both cylindrical and spherical example runs. 

The boundary conditions are defined by the values of "IBTMP" (temperature) and 
"IFREE" (u,v,w components). Both IBTMP and IFREE are arrays of four elements which 
correspond to the four sides of the grid. A heated wall or a no-heated wall are specified by 
setting the respective IBTMP to - 1 or 1 . Similarly, for rigid or free walls, set IFREE to either - 1 
or 1. The values of IFREE and IBTMP do not influence polar conditions for the computations, 
although to obtain good plots they should both be set to 1. 

Printed output and plots 

Selected output is sent to unit 6. The majority of the input stream is written within the 
first 2 pages. After that, diagnostics for the axisymmetric field are given (if not running in 
IRUNCD = 2 mode), and then the waves are written (one per page). Diagnostics include 
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Table A. 1 . GEOSIM Input File Definitions. 


Parameter 

Definition 

NSTEP 

Number of time steps 

ISKIP 

Skip factor for diagnostics (leave zero for default value of 
diagnostics calulated 202 times) 

IRUNCD 

Mode of model integration 

IBTMP 

Temperature boundary conditions 

IFREE 

Momentum boundary conditions 

IPROF 

Profile for initial horizontal temperature gradient 

COND 

Thermal diffusivity 

vise 

Kinematic viscosity 

EXPAN 

Thermal expansivity 

OMEGA 

Rotation rate 

DOMEG 

Used for spin-up cases 

DELT 

Time step 

DTFAC 

Time step factor for temperature 

AXIFAC 

Time step factor for the axisymmetric part 

GAMT 

Optional ADI for temperature (zero for ADI, one otherwise) 

GAMV 

Optional ADI for velocities (same comment) 

RMAX 

Maximum divergence residual (criterion for pressure correction) 

ICONV 

Wave convergence parameter 

ISTART 

Restart option 

PHI1 

Lower meridional boundary (degrees for spherical cases, radians 
for cylindrical) 

PHI2 

Upper meridional boundary (same comment) 
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Table A.2. Memory and Time Requirements for Typical GEOSIM Integrations. 



Symmetric 

only 

Linear 

Wave/Mean 

flow 

Nonlinear 
(1 wave) 

Nonlinear 
(15 waves) 

Option 

1 

2 

4 

5 

5 

30X30 grid 
50X20 grid 

0.39 (62) 
0.32 (65) 

0.44(127) 

0.35(135) 

0.44 (223) 
0.35 (238) 

0.46 (229) 
0.38 (245) 

1.26 (2170) 
1.25 (2429) 


The first number indicates memory requirements in CRAY megawords. The number in 
parenthesis is the number of seconds that were required to run 5000 time steps. 
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integrals of such quantities as kinetic energy, heat fluxes, and energy conversions. The user is 
referred to the subroutines beginning with DIAG... for the definitions of the diagnostics. 

A plot routine (separate from the model) runs on the IBM front-end, is menu driven, and 
plots cross sections of fields (either total or components) in any one of three planes. It also 
allows viewing of certain integrated diagnostics, plotted as a function of time. Because it is 
unlikely that a potential user would be able to transport this code to a different installation, 
documentation is not given here. 
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